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ABSTRACT 


This  report  provides  documentation  and  example  results  for  a computer  model  developed 
to  predict  the  surface  temperature  and  time-of-wetness  (and  time-of-freezing)  of  concrete 
pavements  and  bridge  decks.  The  model  is  based  on  a one-dimensional  finite  difference 
scheme  and  includes  heat  transfer  by  conduction,  convection,  and  radiation.  Environmental 
conditions  are  varied  by  using  the  typical  meteorological  year  datafiles  available  from  the 
National  Renewable  Energy  Laboratory.  Based  on  the  weather  data  available  in  these  files 
and  the  developed  heat  transfer  model,  both  the  surface  temperature  and  the  time-of-wetness 
of  the  concrete  are  computed  for  typical  pavement  and  bridge  deck  structures.  The  time- 
of-wetness  includes  both  precipitation  and  condensation  on  the  concrete  surface.  These 
predictions  will  be  utilized  as  input  into  sorptivity-based  service  life  models  for  sulfate  attack 
and  freeze/thaw  deterioration,  which  is  the  main  goal  of  the  NIST  research  project  being 
funded  by  the  Federal  Highway  Administration. 

Keywords:  Bridge  decks,  building  technology,  computer  modeling,  concrete,  heat  transfer, 
pavements,  surface  temperature,  time-of-wetness. 
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1 Introduction 


As  part  of  a research  project  being  funded  by  the  Federal  Highway  Administration  (FHWA), 
the  National  Institute  of  Standards  and  Technology  (NIST)  is  developing  a computer  pro- 
gram to  predict  the  service  life  of  concrete  pavements  and  bridge  decks  exposed  to  sulfate 
attack  and/or  freeze-thaw  conditions  [1].  For  these  models,  the  major  mechanism  of  trans- 
port of  water  and  sulfate  ions  being  considered  is  sorption  by  the  partially  saturated  concrete. 
The  uptake  of  water  and  the  sulfate  ions  contained  within  it  is  dependent  on  many  variables 
including  the  time-of-wetness  of  the  concrete  surface,  the  saturation  of  the  concrete,  and  its 
surface  temperature.  The  objective  of  this  report  is  to  describe  the  details  of  a computer 
model  for  predicting  these  three  variables  as  a function  of  the  environmental  exposure  for 
hardened  concrete.  Models  for  the  thermal  properties  and  temperature  distribution  of  fresh 
hardening  concrete  have  been  presented  in  the  literature  previously  [2,  3]. 

The  developed  computer  model  (CONCTEMP)  is  based  on  a one-dimensional  finite  dif- 
ference scheme  [4]  and  includes  heat  transfer  due  to  conduction,  convection,  and  radiation. 
The  top  surface  boundary  conditions  (and  the  lower  surface  boundary  conditions  for  bridge 
decks)  are  established  based  on  the  typical  meteorological  years  weather  data  (TMY2DATA) 
provided  by  the  National  Renewable  Energy  Laboratory  (NREL)  [5].  Material  thermal  prop- 
erties and  convection  coefficients  are  taken  from  the  available  literature.  Outputs  from  the 
model  include  the  variation  of  concrete  surface  temperature  during  the  course  of  a year,  a file 
detailing  all  time-of-wetness  events  (start  time,  duration,  concrete  surface  temperature,  and 
exterior  relative  humidity  prior  to  event),  and  a file  detailing  time-of-freezing  events  (start 
time,  duration,  and  minimum  surface  temperature  achieved  during  freezing).  These  files  have 
been  created  for  twelve  representative  geographical  locations  throughout  the  United  States. 
A complete  listing  of  the  computer  code  used  to  estimate  the  concrete  surface  temperature 
for  pavement  structures  can  be  found  in  Appendix  A.  The  program  has  been  written  in  the  C 
programming  language.  It  has  been  written  in  a modular  fashion  with  separate  subroutines 
for  reading  in  the  input  material  properties  and  for  updating  the  temperatures  of  all  nodes 
within  the  1-D  structure.  The  log  files  for  wetting  and  freezing  events  are  created  within  the 
main  program  routine,  which  is  also  responsible  for  accumulating  time  throughout  the  year 
and  reading  in  the  hourly  weather  data  values. 

DISCLAIMER 

This  software  was  developed  at  the  National  Institute  of  Standards  and  Technology  by 
employees  of  the  Federal  Government  in  the  coursse  of  their  official  duties.  Pursuant  to  title 
17  Section  105  of  the  United  States  Code  this  software  is  not  subject  to  copyright  protection 
and  is  in  the  public  domain.  CONCTEMP  is  an  experimental  system.  NIST  assumes  no 
responsibility  whatsoever  for  its  use  by  other  parties,  and  makes  no  guarantees,  expressed 
or  implied,  about  its  quality,  reliability,  or  any  other  characteristic.  We  would  appreciate 
acknowledgement  if  the  software  is  used. 

The  U.S.  Department  of  Commerce  makes  no  warranty,  expressed  or  implied,  to  users 
of  CONCTEMP  and  associated  computer  programs,  and  accepts  no  responsibility  for  its 
use.  Users  of  CONCTEMP  assume  sole  responsibility  under  Federal  law  for  determining  the 
appropriateness  of  its  use  in  any  particular  application;  for  any  conclusions  drawn  from  the 
results  of  its  use;  and  for  any  actions  taken  or  not  taken  as  a result  of  analyses  performed 
using  these  tools. 
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Users  are  warned  that  CONCTEMP  is  intended  for  use  only  by  those  competent  in  the 
field  of  cement-based  materials  and  is  intended  only  to  supplement  the  informed  judgment 
of  the  qualified  user.  The  software  package  is  a computer  model  which  may  or  may  not  have 
predictive  value  when  applied  to  a specific  set  of  factual  circumstances.  Lack  of  accurate  pre- 
dictions by  the  model  could  lead  to  erroreous  conclusions  with  regard  to  materials  selection 
and  design.  All  results  should  be  evaluated  by  an  informed  user. 

INTENT  AND  USE 

The  algorithms,  procedures,  and  computer  programs  described  in  this  report  constitute 
a methodology  for  predicting  the  surface  temperature  and  time-of-wetness  of  concrete  un- 
der typical  weathering  conditions.  They  have  been  compiled  from  the  best  knowledge  and 
understanding  currently  available,  but  have  important  limitations  that  must  be  understood 
and  considered  by  the  user.  The  program  is  intended  for  use  by  persons  competent  in  the 
field  of  cement-based  materials  and  with  some  familiarity  with  computers.  It  is  intended  as 
an  aid  in  the  materials  selection,  optimization,  and  design  process. 


2 Heat  Transfer  Model 


The  basic  geometrical  configurations  for  concrete  pavements  and  bridge  decks  considered  in 
this  study  are  shown  in  Figure  1.  For  a pavement,  the  concrete  thickness  is  considered  to 
be  0.2  m,  with  a 0.2  m thick  layer  of  soil  (subbase)  immediately  beneath  it.  For  the  bridge 
deck,  both  surfaces  of  the  concrete  are  assumed  to  be  exposed  to  the  environment,  which  is 
generally  the  case  when  temporary  wooden  forms  are  used  [6,  7].  While  often  the  bottom 
surface  is  protected  from  the  environment  by  steel  forms  which  remain  in  place  [6,  7],  the 
thickness  of  these  forms  (20  gauge  or  0.9  mm)  is  such  that  their  contribution  to  heat  transfer 
should  be  minimal.  For  the  bridge  decks,  it  is  further  assumed  that  the  convection  coefficient 
for  heat  transfer  is  the  same  for  both  the  top  and  bottom  surfaces,  and  that  no  heat  transfer 
by  radiation  (incoming  sunlight  or  emitted  radiation)  occurs  at  the  bottom  surface.  For  the 
pavements,  it  is  further  assumed  that  at  a depth  of  0.2  m,  the  soil  temperature  is  constant 
at  a value  of  13  °C.  (Variation  of  this  soil  temperature  did  not  significantly  influence  the 
model  results  in  this  study). 

Both  systems  shown  in  Figure  1 are  modelled  using  a one-dimensional  finite  difference 
grid  with  a spatial  resolution  of  20  mm  (e.g.,  ten  nodes  within  the  concrete).  At  the  top 
surface,  four  modes  of  heat  transfer  are  considered:  conduction  into  the  concrete,  convection, 
solar  absorption,  and  grey-body  irradiation  to  the  surroundings.  Specifically,  the  heat  flow 
contribution  (in  ^)  due  to  conduction  is  given  by: 

Qcond,  — kconc  * 

where  kconc  is  the  thermal  conductivity  of  the  concrete  in  m^,  To  and  T\  are  the  surface 
temperature  and  internal  temperature  at  the  first  internal  node,  respectively,  and  Ax  is  the 
node  spacing  (20  mm).  For  convection,  the  heat  flow  is  given  by: 


(Tp  - Ti) 

Ax 


a) 


Qconv  h-conv  ^ (To  Tam^en^) 


(2) 


where  Tam^ent  is  the  ambient  temperature  and  hconv  is  the  convection  coefficient  in 
While  several  empirical  relationships  to  estimate  convection  coefficients  are  available  in  the 
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Figure  1:  Basic  configurations  of  one-dimensional  heat  transfer  model  for  concrete  pavements 
and  bridge  decks. 
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literature  [2],  for  this  study,  the  convection  coefficient  was  calculated  based  on  the  wind 
speed  available  in  the  weather  data  files  and  the  following  equations  used  in  the  commercially 
available  FEMMASSE  system  [8]: 

hconv  = 5.  6 + 4.0vwind  for  vwind  < 5 m/s  (3) 

hconv  = 7.2  x for  vwind  > 5 m/s 

where  vWind  is  the  measured  wind  speed  in  m/s. 

For  radiative  heat  transfer  at  the  top  concrete  surface,  two  contributions  are  considered. 
The  first  is  the  radiation  absorbed  from  the  incoming  sunlight.  A simple  equation  for  the 
incoming  heat  flow  due  to  this  source  is  given  by  [2]: 

Qsun  = 'Jabs  * Qinc  (4) 

where  Qinc  is  the  incident  solar  radiation  (W/m2)  and  7^  is  the  solar  absorptivity  of  the 
concrete.  For  this  study,  the  incident  solar  radiation  was  taken  directly  from  the  weather 
data  files  [5],  while  a value  of  0.65  was  used  for  the  solar  absorptivity  [8].  Second,  one  must 


3 


consider  the  emission  of  radiation  from  the  “warm”  concrete  to  the  night  sky.  This  heat  flow 
can  be  estimated  as: 

Qsky  ==  ere  x (Tok  — Tsky)  (5) 

where  a is  the  Stefan-Boltzmann  constant  (5.669  x 10-8  W/(m2  °C4)),  e is  the  emissivity 
of  the  concrete,  TQK  is  the  concrete  surface  temperature  (in  K),  and  Tsky  is  the  calculated 
sky  temperature  also  in  K.  For  this  study,  a value  of  0.9  was  used  for  the  concrete  emissivity 
[2,  8,  9].  The  sky  temperature  was  estimated  based  on  an  algorithm  presented  by  Walton 
[10]  using  the  following  series  of  equations: 

Tsky  — es  4 X T ambient  {F  in  K ) (6) 

where  the  sky  emissivity,  es  is  given  by: 

es  = 0.787  + 0.764  x * Fclmd  (7) 

where  Tdew  is  the  dewpoint  temperature  in  K and  with  the  cloud  cover  factor,  Fd^d,  as: 

Fdoud  = 1.0  + 0.0247V  - 0.00357V2  + 0.000287V3  (8) 

where  TV  is  the  “tenths  cloud  cover”,  taking  values  between  0.0  and  1.0.  Because  the  above 
equations  for  surface  irradiation  are  significantly  different  from  approximations  used  else- 
where [2],  it  is  worth  noting  that  the  equations  provided  here  have  been  successfully  used  in 
a variety  of  predictive  models  of  relevance  to  the  construction  community  [10,  11]. 

The  above  equations  are  employed  in  a general  finite  difference  solution  to  one-dimensional 
heat  transfer  [4].  The  time  step  in  the  finite  difference  scheme  is  established  based  on  the 
layer  thickness  (20  mm)  and  the  thermal  properties  of  the  materials,  to  ensure  numerical 
convergence  of  the  solution  [4].  The  ambient  environmental  conditions  at  any  specific  time 
are  calculated  by  linear  interpolation  of  the  hourly  values  available  from  the  weather  files. 

2.1  Material  Properties 

The  additional  material  properties  needed  for  the  heat  transfer  model  are  provided  in  Table  1. 
The  values  for  the  soil  were  taken  directly  from  the  Final  Report  for  the  FHWA  HIPERPAV 
project  [2].  For  concrete,  a variety  of  values  for  the  properties  are  available  in  the  literature 
[2,  3,  4,  12]  and  the  values  provided  in  Table  1 are  considered  to  be  reasonable  “mean” 
estimates.  Of  course,  it  is  well  recognized  that  each  of  these  material  properties  will  vary 
with  the  specific  concrete  under  consideration.  If  specific  values  were  measured  for  a concrete 
of  interest,  the  heat  transfer  model  could  be  easily  executed  for  these  values,  since  they  are 
typically  provided  in  an  input  file  to  the  main  program. 

2.2  Environmental  Conditions 

When  considering  the  service  life  of  a field  concrete,  the  environmental  conditions  to  which 
the  concrete  is  exposed  are  equally  as  important  as  its  transport  and  mechanical  properties. 
Naturally,  these  conditions  are  highly  variable  thoughout  the  United  States.  For  the  current 
study,  realistic  environmental  exposures  were  obtained  by  using  the  typical  meteorological 
year  weather  data  files  available  from  NREL  [5].  These  data  files  provide  typical  hourly 
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Table  1:  Material  Properties  for  Heat  Transfer  Model 


Heat  Capacity 

Thermal  Conductivity 

Density 

Material 

Cp  (ko  °c) 

* (mwJ 

p&) 

Concrete 

1000 

1.5 

2350 

Soil 

800 

0.3 

1600 

weather  data  for  over  200  geographical  locations  within  the  United  States,  each  geographical 
location  being  characterized  by  a five  digit  numerical  code.  For  the  heat  transfer  model 
described  above,  the  following  TMY2DATA  parameters  are  utilized:  ambient  temperature 
(' Tambient ),  dewpoint  temperature  (7detw),  ambient  relative  humidity,  wind  speed  ( rwd ),  pre- 
cipitation, cloud  cover  (iV),  and  incident  global  horizontal  radiation  (Qinc)- 

As  the  heat  transfer  model  executes,  two  types  of  events  are  logged,  namely  wetting  and 
freezing.  A wetting  event  may  be  due  either  to  a precipitation  event  (rain,  drizzle,  sleet,  or 
hail)  in  the  weather  data  file  or  condensation  when  the  computed  surface  temperature  of  the 
concrete  falls  below  the  current  dewpoint  temperature  (but  is  still  above  0 °C).  It  should  be 
noted  that  snow  is  not  considered  a wetting  event  as  we  are  assuming  that  the  snow  will  be 
successfully  removed  from  the  pavement  or  bridge  deck  before  substantial  melting/wetting 
occurs.  A freezing  event  is  defined  by  the  concrete  surface  temperature  falling  below  0 °C. 
This  value  was  chosen  as  a conservative  estimate,  since  the  freezing  point  of  pore  water  in 
concrete  is  typically  depressed  by  the  presence  of  ionic  species  and  the  small  sizes  of  the 
capillary  and  gel  pores  in  the  concrete.  When  a wetting  event  occurs,  the  following  values 
are  logged  to  a file:  the  initiation  time  (hours  since  the  beginning  of  the  year),  the  duration 
of  the  event  in  hours,  the  average  of  the  concrete  surface  temperatures  at  the  initiation 
and  the  termination  of  the  event,  and  the  ambient  relative  humidity  at  the  initiation  of  the 
event.  For  a freezing  event,  only  the  initiation  time,  duration,  and  minimum  achieved  surface 
temperature  are  logged  to  a freezing  event  data  file. 


3 Example  Results 

3.1  Concrete  Surface  Temperature  vs.  Time  of  Year 

Three  examples  of  predicted  concrete  surface  temperatures  as  a function  of  geographical 
location  and  season  of  the  year  are  provided  in  Figures  2,  3,  and  4.  Little  data  is  available  in 
the  literature  to  validate  these  predictions,  but  the  data  of  both  Andrade  [13]  and  Basheer 
and  Nolan  [14]  support  the  general  trends  shown  in  the  three  figures.  Specifically,  during  the 
day,  the  concrete  surface  temperature  generally  rises  above  the  ambient  temperature  due  to 
the  incoming  solar  radiation.  At  night,  the  concrete  temperature  falls  due  to  radiation  from 
the  concrete  surface  to  the  sky,  sometimes  falling  below  the  ambient  air  temperature  and 
occasionally  falling  below  the  dewpoint.  This  effect  is  more  significant  for  bridge  decks  than 
for  pavements  since  there  is  no  soil  sublayer  providing  a thermal  mass  beneath  the  concrete 
in  the  case  of  bridge  decks. 
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Figure  2: 


Temperature  predictions  for  a concrete  pavement  surface 


in  Seattle,  Washington. 
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Figure  3:  Temperature  predictions  for  a concrete  pavement  surface  in  Tucson,  Arizona. 
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Figure  4:  Temperature  predictions  for  a concrete  bridge  deck  surface  in  Alpena,  Michigan. 
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3.2  Time-of- Wetness  and  Time-of- Freezing  Predictions 

The  wetting  and  freezing  events  for  12  representative  geographical  locations  are  summarized 
in  Tables  2 and  3.  It  is  easily  observed  that  there  is  considerable  variability  in  time-of-wetness 
and  number  of  freeze-thaw  cycles  across  the  United  States.  Of  the  12  locations  investigated 
to  date,  Tucson,  AZ  is  the  driest  with  the  concrete  pavement  being  wet  only  128  h per  year. 
Conversely,  Seattle,  WA  is  the  wettest  climate  with  the  pavement  being  wet  1698  h per 
year.  Obviously,  the  time-of-wetness  and  number  of  freeze-thaw  cycles  will  have  a significant 
influence  on  the  durability  of  concrete  pavements  and  bridge  decks  in  any  given  geographical 
location.  This  influence  will  be  further  developed  in  the  sorptivity-based  service  life  models 
for  sulfate  attack  and  freeze-thaw  deterioration  [1]  currently  being  developed  as  one  of  the 
final  stages  of  this  research  project. 

For  freezing  events,  there  are  always  more  freeze-thaw  cycles  for  the  bridge  decks  than  for 
the  pavements  (in  agreement  with  the  ubiquitous  posting  of  road  signs  indicating  that  “bridge 
freezes  before  roadway”).  Not  surprisingly,  Tampa,  FL  experiences  the  fewest  freeze-thaw 
cycles,  with  only  the  bridge  decks  experiencing  any  cycles  in  a representative  year.  Cheyenne, 
WY  experiences  the  most  freeze-thaw  cycles,  with  over  125  cycles  for  both  the  pavements 
and  bridge  decks  being  observed  in  a representative  year.  In  comparing  bridge  decks  to 
pavements,  Baltimore,  MD  exhibits  the  greatest  difference,  with  the  bridge  deck  concrete 
experiencing  21  more  freeze-thaw  cycles  per  year  than  a comparable  pavement.  This  is  most 
likely  due  to  the  somewhat  moderate  winter  climate  in  Baltimore,  with  the  ambient  often 
near  0 °C,  so  that  the  bridge  decks  often  experience  freezing  while  the  nearby  pavements  are 
just  above  the  freezing  point. 


Table  2:  Wetting  and  Freezing  Events  for  Concrete  Pavements 


U.S.  City 

NREL  code 

Number  of 
wetting  events 

Total  h/yr 
wet 

Average  RH 
before  wetting 

Number  of 
F-T  cycles 

Kansas  City,  MO 

03947 

173 

448 

83.7 

79 

Tampa,  FL 

12842 

244 

613 

86.0 

0 

Lubbock,  TX 

23042 

138 

348 

83.6 

61 

Tucson,  AZ 

23160 

78 

128 

66.3 

9 

Cheyenne,  WY 

24018 

117 

224 

71.5 

126 

Pierre,  SD 

24025 

117 

190 

85.3 

92 

Seattle,  WA 

24233 

407 

1698 

81.6 

25 

Fresno,  CA 

93193 

139 

393 

86.2 

14 

Baltimore,  MD 

93721 

238 

948 

84.7 

83 

Bridgeport,  CT 

94702 

229 

835 

87.2 

90 

Alpena,  MI 

94849 

218 

680 

86.1 

102 

Waterloo,  IA 

94910 

209 

551 

86.7 

72 
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Table  3:  Wetting  and  Freezing  Events  for  Concrete  Bridge  Decks 


U.S.  City 

NREL  code 

Number  of 
wetting  events 

Total  h/yr 
wet 

Average  RH 
before  wetting 

Number  of 
F-T  cycles 

Kansas  City,  MO 

03947 

179 

490 

84.1 

81 

Tampa,  FL 

12842 

268 

639 

86.8 

4 

Lubbock,  TX 

23042 

148 

413 

84.3 

71 

Tucson,  AZ 

23160 

78 

141 

66.3 

16 

Cheyenne,  WY 

24018 

125 

252 

72.9 

131 

Pierre,  SD 

24025 

123 

229 

85.7 

100 

Seattle,  WA 

24233 

413 

1812 

81.7 

34 

Fresno,  CA 

93193 

161 

521 

87.0 

20 

Baltimore,  MD 

93721 

254 

972 

85.3 

104 

Bridgeport,  CT 

94702 

229 

860 

87.0 

104 

Alpena,  MI 

94849 

239 

758 

86.8 

107 

Waterloo,  IA 

94910 

212 

623 

86.8 

86 

4 Summary 

A simple  one-dimensional  heat  transfer  model  has  been  applied  to  estimating  the  surface 
temperature  of  concrete  pavements  and  bridge  decks.  Using  representative  weather  data 
from  the  NREL,  the  results  indicate  considerable  geographical  variability  not  only  in  the 
concrete  surface  temperatures,  but  also  in  the  time-of-wetness  and  number  of  freeze-thaw 
cycles.  This  indicates  that  realistic  estimates  of  service  life  can  only  be  made  when  detailed 
consideration  is  given  to  the  environment  specific  to  the  field  concrete  in  question.  The  results 
of  this  model  will  be  utilized  as  input  to  a Windows-based  computer  program  for  sorptivity- 
based  service  life  predictions  for  the  cases  of  sulfate  attack  and  freeze-thaw  deterioration, 
currently  under  development  at  NIST. 
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Appendix 


A Source  Code  Listing  for  Heat  Transfer  Model  for 
Concrete  Pavements 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <math.h> 

#define  STEFAN  5.669e-8  /*  Stef an-Boltzmann  constant  */ 

#define  EMISS  0.90  /*  Emissivity  of  concrete  surface  */ 

#define  SOLARABS  0.65  /*  Solar  absorptivity  of  concrete  */ 

/*  */ 

/*  Program:  Tempcalc.c  */ 

/*  Program  to  calculate  concrete  surface  temperature  during  a */ 

/*  yearly  weather  cycle  (based  on  a 1-D  finite  difference  solution  */ 

/*  */ 

/*  Inputs:  File  inputs. par  with  material  properties  of  concrete  */ 

/*  and  soil.  */ 

/*  Weather  data  file  XXXXX.tm2  from  NREL  database  */ 

/*  Outputs:  Files  with  */ 

/*  1)  surface  temperature  history  (temphist .out)  */ 

/*  2)  time-of-wetness  history  (wettime. out)  */ 

/*  3)  time-of -freezing  history  (freezetime . out)  */ 

/*  */ 

/*  Programmer:  Dale  P.  Bentz  */ 

/*  Last  revision:  July  19,  2000  */ 

/*  Effort  funded  by  the  Federal  Highway  Administration  */ 

/*  */ 

/*******************  + ****  + + **4:**  + + + ***=(C5|c*  + 4:*******************=t:**********/ 

/*  Global  variables  */ 

double  tempnew [22] , tempold [22] , t _out=0 . 0 , sum_t , tk_sky ; 
double  conc_heat_cap,conc_thermal_k,conc_density ; 
double  h_coef f =0 . 0 , alpha_conc , t _cum=0 . 0 ; 

double  soil_heat_cap,soil_thermal_k, soil_density ,alpha_soil; 

/*  Routine  to  read  in  input  property  parameters  for  concrete  and  soil  */ 

/*  Called  by:  main  program  */ 

/*  Calls:  none  */ 
read_inputs(){ 

FILE  *parfile; 
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parfile=f open ("inputs .par" ,"r") ; 

/*  Note  that  all  inputs  are  in  SI  units  */ 
f scanf  (parf  ile,  "‘/.If  " ,&conc_heat_cap)  ; 

printf ("Concrete  heat  capacity  is  %lf  J/kg/C\n" ,conc_heat_cap) ; 
f scanf (parf ile, "%lf " ,&conc_thermal_k) ; 

printf ("Concrete  thermal  conductivity  is  %lf  W/m/C\n" ,conc_thermal_k) ; 
f scanf (parf ile, "%lf " ,&conc_density) ; 

printf ("Concrete  density  is  %lf  kg/m/m/m\n" , conc_density) ; 
alpha_conc=conc_thermal_k/conc_heat_cap/conc_density; 
printf  ("Calculated  concrete  alpha  value  is  7,.121f  m*m/s\n" ,alpha_conc) ; 
f scanf  (parf  ile,  "*/elf  " ,&soil_heat_cap)  ; 

printf  ("Soil  heat  capacity  is  7#lf  J/kg/C\n" ,soil_heat_cap) ; 
f scanf  (parf  ile,  "’/.If  " ,&soil_thermal_k) ; 

printf  ("Soil  thermal  conductivity  is  '/.If  W/m/C\n" ,soil_thermal_k) ; 

f scanf  (parf  ile,  "Velf"  ,&soil_density)  ; 

printf  ("Soil  density  is  7»lf  kg/m/m/m\n"  ,soil_density) ; 

alpha_soil=soil_thermal_k/soil_heat_cap/soil_density ; 

printf  ("Calculated  soil  alpha  value  is  °/,.121f  m*m/s\n" ,alpha_soil) ; 

f close (parf ile) ; 

> 

/*  Routine  to  update  all  temperatures  in  the  finite  difference  grid  */ 

/*  Called  by:  main  program  */ 

/*  Calls:  none  */ 

upd_temp(q_in , t_amb , t_dew , cloud_cover , delta.x , delta.t) 

double  q_in , t_amb , t_dew , cloud.cover , delta_x , delta_t ; 

/*  q_in  is  in  W/m*m 

t_amb  and  t_dew  are  in  degrees  Celsius 
cloud_cover  is  in  tenths  (0. 1-1.0) 
delta_x  is  in  meters 
delta_t  is  in  seconds  */ 

{ 

int  i ; 

FILE  *tempf ile ; 

double  q_sun , q_sky , q_conv , q_cond , q_st or , ql , q2 ; 
double  tk_amb,tk_dew,e_s, cloud_fact ; 

tk_amb=t_amb+273. 15;  /*  Conversion  to  Kelvin  */ 
tk_dew=t_dew+273. 15;  /*  Conversion  to  Kelvin  */ 

/*  Sky  temperature  computed  based  on  algorithm  of  Walton  et  al.  (TARP)  */ 
cloud_f act=l . 0+0 . 024*cloud_cover-0 . 0035*cloud_cover*cloud_cover ; 
cloud_f act+=0 . 00028*cloud_cover*cloud_cover*cloud_cover ; 

/*  Compute  the  sky  emissivity  e_s  */ 

e_s= (0 . 787+0 . 764*log(tk_dew/273 . ) ) *cloud_f act ; 

/*  Now  estimate  the  sky  temperature  tk_sky  */ 
tk_sky-sqrt (sqrt (e_s) ) *tk_amb ; 
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/*  Compute  the  four  heat  flow  contributions  at  the  top  surface  */ 
q_sun=q_in*SOLARABS ; 

q_sky=STEFAN*EMISS* (pow ( (t empold [0] +273 . 15) ,4.0) -pow (tk_sky ,4.0)); 
q_conv=h_coeff * (tempold[0] -t_amb) ; 

q_cond= (conc_thermal_k/delta_x) * (tempold [0] -tempold [1] ) ; 
q_stor=q_sun-q_sky-q_conv-q_cond; 

/*  Update  the  surface  temperature  */ 

tempnew [0] =t empold [0] +q_stor*2 . *delta_t/delta_x/conc_density/conc_heat_cap ; 
/*  First,  the  top  concrete  layer  */ 

/*  Equations  based  on  those  presented  in  Heat  Transfer  by  Holman  */ 
f or(i=l  ;i<10;  i++M 

tempnew [i] = ( (alpha_conc*delta_t) / (delta_x*delta_x) ) * 

(tempold [i+l]+tempold[i-l] )+(l .0-2. *alpha_conc*delta_t/ 
(delta_x*delta_x))*tempold[i] ; 

> 

/*  Second,  the  interface  between  concrete  and  soil  */ 
ql= (conc_thermal_k/delta_x) * (tempold [10] -tempold [9] ) ; 
q2= (soil_thermal_k/delta_x) * (tempold [10] -tempold [11] ) ; 
tempnew [10] =tempold [10] - (ql+q2) *2 . *delta_t/delta_x/ 
(conc_density*conc_heat_cap+soil_density*soil_heat_cap) ; 

/*  Last,  the  soil  subbase  */ 
f or  (i=ll;  i<20 ; i++M 

tempnew  [i]  = ( (alpha_soil*delta_t) / (delta_x*delta_x) ) * 

(tempold [i+1] +tempold [i-1] ) + ( 1 . 0-2 . *alpha_soil*delta_t/ 
(delta_x*delta_x))*tempold[i] ; 

> 

/*  Last  soil  node  is  constant  at  a temperature  of  13  C */ 
tempnew [20] =13 . 0 ; 

/*  Copy  the  new  temperatures  over  the  old  ones  to  prepare  for  the  next  update  */ 
f or (i=0 ; i<=20 ; i++) { 

tempold  [i]=tempnew[i]  ; 

} 

/*  Output  temperature  results  once  per  hour  */ 
if ( (t _cum-t_out ) >=3600 . ) { 

tempf ile=f open ("temphist .out" ,"a") ; 
f printf (tempf ile , "% . 21f  */,.21f  %.21f  %.21f  ", 

t_cum/3600. ,t_amb,t_dew,tk_sky-273. 15) ; 
t_out=t_cum; 

/*  Output  surface  temperature  */ 
f or (i=0 ; i<=l ; i++) { 

f printf  (tempf  ile,"  V.lf  ",  tempnew  [i] ) ; 

} 
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/*  And  temperature  for  every  fifth  node  into  concrete  and  soil  */ 
f or (i=5 ; i<=20 ; i+=5) { 

f printf (tempf ile , "%lf  " , tempnew [i] ) ; 

> 

f printf (tempf ile, "\n") ; 
f close (tempf ile) ; 

} 

> 

/*  Main  program  to  read  in  weather  data  file  and  call  routines  to  */ 

/*  set  material  properties  and  update  temperatures  */ 

/*  Calls:  read_inputs  and  upd_temp  */ 
main(){ 

FILE  *infile,*wetfile,*frfile; 

char  wba[80]  , city [80]  .state [80]  ,latd[10]  ,longd[10]  ; 
char  f ilein[80] ; 
register  int  i,j; 

int  tzone , lat 1 , lat2 , longl , long2 , elev ; 
int  wetf lag=0 , duration ,numcum=0 , f reezef lag=0 ; 
long  int  timeorg.timef in; 
long  int  timef org.timeff in; 

int  year .month , day , hour , exhr , exdnr , glhr , glhrf lu , dnr , dnrf lu ; 

int  df hr , df hrf lu , glhi , glhif lu , dhi , dhif lu , dni , dnif lu ; 

int  zi , zif lu , t skycov , t skycovf lu , opqskycov , opqskycovf lu , tempdry , tdf lu ; 

int  tempwet , twf lu , rh , rhf lu , atmp , atmpf lu , wdir , wdirf lu , wspd , wspdf lu ; 

int  vis,visflu,ceilflu,days[13] ; 

long  int  ceil.durcum; 

double  qsun , temp_amb , temp_dew , del_x , del_t , del_t_soil , tf act ; 
double  temp_amb_old , wspd_old , temp_dew_old , wspd_new ; 
double  temp_amb_cur , wspd.cur , temp_dew_cur , tf min=0 . 0 ; 
float  rhsum=0 . 0 , rhave , rhinit , twinit=0 . 0 ; 

int  pwO , pwl , pw2 , pw3 , pw4 , pw5 , pw6 , pw7 , pw8 , pw9 , prewat , pr ewatf lu ; 
char  glhrf Is , dnrf Is , df hrf Is , glhif Is , dhif Is , dnif Is , rest [240] ; 
char  zif Is , t skycovf Is , opqskycovf Is , tdf Is , twf Is , rhf Is , atmpf Is .wdirf Is ; 
char  wspdf Is , visf Is , ceilf Is .prewatf Is ; 

read_inputs() ; 

del_x=0.02;  /*  2 cm  spacing  in  nodes  for  40  cm  total  */ 

durcum=0 ; 

/*  Compute  cumulative  days  per  month  */ 

days[l]=0; 

days [2] =31; 

days [3] =days [2] +29 ; 

days  [4]  =days  [3]  +31 ; 

days  [5]  =days  [4]  +30 ; 

days  [6]  =days  [5]  +31 ; 
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days [7] =day s [6] +30 ; 
days  [8]  =days  [7]  +31 ; 
days [9] =days [8] +31 ; 
days [10] -days [9] +30 ; 
days [11] =days [10]  +31 ; 
days [12] =days [11] +30 ; 

/*  Get  name  of  weather  data  file  to  use  for  this  simulation  */ 

printf ("Enter  name  of  file  to  read  \n"); 

scanf  ("7oS"  ,filein)  ; 

printf  ("'/»s\n"  ,filein)  ; 

inf ile=fopen(filein, "r") ; 

/*  Open  log  files  for  recording  time-of -wetness  and  time-of -freezing  */ 
wetf ile=f open ( "wettime . out " , " w" ) ; 
frf ile=fopen("freezetime .out" , "w") ; 

/*  Read  in  the  header  for  this  weather  data  file  and  echo  values  */ 
f scanf  (infile,  "7.s  7,s  Is  7.d  Is  7.d  7,d  Is  %d  7.d  7.d", 

wba , city , state , fetzone , latd , felat 1 , &lat 2 , longd , felongl , &long2 , feelev) ; 
printf  ("Number  is  7oS\n",wba); 
printf  ("City  is  7.s\n"  ,city) ; 
printf  ("State  is  7«s\n" , state)  ; 
printf  ("Time  zone  is  7od\n"  ,tzone)  ; 

/*  Set  time  step  in  finite  difference  routine  */ 

/*  based  on  thermal  properties  of  concrete  and  soil  */ 
del_t=del_x*del_x/4 . /alpha_conc ; 
del_t_soil=del_x*del_x/4 . /alpha_soil ; 
if (del_t_soil<del_t) {del_t=del_t_soil ; > 

/*  Scale  time  step  to  a factor  of  3600  (seconds  in  an  hour)  */ 
if (del_t<50.0){ 

del_t=10 . ; 

} 

else  if ( (del_t>50 . 0) && (del_t<100 . 0) ) { 
del_t=50. 0; 

> 

else  if  ( (del_t>100 . 0)&&  (del_t<200 . 0) ) ■[ 
del_t=100.0; 

> 

else  if ( (del_t>200 . 0)&& (del_t<300 . 0) ) { 
del_t=200.0; 

} 

else  if ( (del_t>300.0)&&(del_t<400.0)){ 
del_t=300.0; 

} 

else  if (del_t>400.0){ 
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del_t=400.0; 


> 

printf ("del_t  is  '/.If  \n",del_t); 
fflush(stdout) ; 

/*  Read  in  the  first  record  */ 

/*  Format  statement  is  directly  from  TMY2DATA  (NREL)  Manual  */ 


/oio/ol  U/e*IU/o  lD/«  iU/otU/j  X O /oJ.U/o^U/0  iO/ol  U/oXU/o  X O /oXU/0TU/e  x Q /olU/olU/o  X O /o  x U/eJU/o  X O /0  lU/o^U/ol  O /o  XU/o 

d,/.ls0/.ld'/.3d'/.ls'/.ld'/,4d0/,ls'/.ld,/»51d'/.ls’/»ld7.1d7»ld'/.ldy.ld7.1d0/.ld,/,ld,/.ld,/old'/»ld’/,3d’/,ls7.1d7.s" , 
feyear , femonth , feday , fehour , feexhr , feexdnr , feglhr , feglhrf Is , feglhrf lu , fednr , fednrf Is , 
fednrf lu , fedf  hr , fedf hrf Is , &df hrf lu , feglhi , feglhif Is , feglhif lu , &dni , fednif Is , fednif lu , 
&dhi , fedhif Is , fedhif lu, fezi , &zif Is , fezif lu, fetskycov , fetskycovf Is , fetskycovf lu , 
feopqsky cov , feopqskycovf Is , feopqskycovf lu , &t empdry , fetdf Is , &tdf lu , &t empwet , 
fetwf Is , fetwf lu , &rh , ferhf Is , ferhf lu ,&atmp , featmpf Is , featmpf lu , &wdir , fewdirf Is , 
fewdirf lu , fewspd , fewspdf Is , fewspdf lu , fevis , fevisf Is , fevisf lu , &ceil , &ce ilf Is , 

&ce ilf lu , &pwO , &pwl , &pw2 , &pw3 , &pw4 , &pw5 , &pw6 , &pw7 , &pw8 , &pw9 , fcprewat , 
&prewatfls,&prewatf lu,rest) ; 

/*  Extract  the  relevant  values  */ 
t emp_amb_old= (double) tempdry/ 10 . ; 
wspd_old= (double) wspd/10 . ; 
temp_dew_old= (double) tempwet/10. ; 

/*  Initialize  temperature  data  to  ambient  temperature  value  */ 
f or ( j =0 ; j <=20 ; j ++) { 

tempold [ j ] =temp_amb_old ; 

> 

tempold [20] =13 . 0 ; 


/*  Now  process  each  weather  record  in  turn  */ 
f or (i=l ; i<8760 ; i++) { 


feyear , femonth , feday , fehour , feexhr , feexdnr , feglhr , feglhrf Is , feglhrf lu , fednr , fednrf Is , 
fednrf lu , fedf hr , fedf hrf Is , fedf hrf lu , feglhi , feglhif Is , feglhif lu , fedni , fednif Is , 
fednif lu , fedhi , fedhif Is , fedhif lu , fezi , fezif Is , fezif lu , fetskycov , fetskycovf Is , 
fetskycovf lu , feopqskycov , feopqskycovf Is , feopqskycovf lu , &t empdry , fetdf Is , fetdf lu , 
fetempwet , fetwf Is , fetwf lu , ferh , ferhf Is , ferhf lu , featmp , featmpf Is , featmpf lu , fewdir , 
fewdirf Is , fewdirf lu , fewspd , fewspdf Is , fewspdf lu , fevis , fevisf Is , fevisf lu , feceil , 
fece ilf Is , feceilf lu , fepwO , fepwl , &pw2 , &pw3 , &pw4 , &pw5 , &pw6 , &pw7 , &pw8 , &pw9 , feprewat , 
&prewatfls,feprewatflu,rest) ; 

/*  Extract  the  relevant  values  */ 
temp_amb= (double) tempdry/10. ; 
temp_dew= (double) tempwet/10. ; 
wspd_new= (double) wspd/10. ; 
qsun=glhr ; 
sum_t=0.0; 
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/*  Now  process  one  hours  worth  of  data  */ 
while (sum_t <3600. ){ 

/*  Update  the  hourly  and  cumulative  times  */ 
sum_t+=del_t ; 
t_cum+=del_t ; 
tfact=(sum_t)/3600. ; 

/*  Linearly  interpolate  the  temperatures  and  wind  speeds  */ 

/*  within  any  1 h period  */ 

temp_amb_cur=temp_amb_old+tf act* (temp_amb-temp_amb_old) ; 
temp_dew_cur=temp_dew_old+tf act* (temp_dew-temp_dew_old) ; 
wspd_cur=wspd_old+tf act* (wspd_new-wspd_old) ; 

/*  Convection  coefficient  based  on  wind  speed  from  FEMMASSE  */ 

/*  of  Intron  in  The  Netherlands  (3/00)  */ 
if (wspd_cur<5 . 0) { 

h_coef f =5 . 6+4 . 0*wspd_cur ; 

> 

else{ 

h_coef f =7 . 2*pow (wspd_cur ,0.78); 

> 

/*  Call  the  routine  to  update  all  temperatures  */ 
upd_temp(qsun,temp_amb_cur,temp_dew_cur , (double) tskycov/10. ,del_x,del_t) ; 
> 

/*  Now  look  for  freezing  and  wetting  events  */ 
temp_amb_old=temp_amb ; 
temp_dew_old=temp_dew ; 
wspd_old=wspd_new ; 

/*  Freeze  events  */ 

if ( (f reezef lag==0)&& (tempold [0] < (0 . 0) ) ) { 
freezef lag=l; 

timef org=(day*24+days [month] *24+hour) ; 
tfmin=tempold[0] ; 

> 

/*  Log  minimum  temperature  achieved  during  freezing  */ 
if ( (freezef lag==l)&&(tf min>tempo Id [0] )){tfmin=tempold[0] ;} 
if ( (f reezef lag==l)&& (tempold [0]  >= (0 . 0) ) ) { 
freezeflag=0; 

timef f in=(day*24+days [month] *24+hour) ; 
dur at i on=t imef f in-t ime f o rg ; 

fprintf  (frfile,"*/#ld  */»d  '/.If \n" , timef org, duration, tfmin)  ; 

> 

/*  Wet  if  rain  (precipitation)  or  dew  (condensation)  event  */ 
if ( ( (wetf lag==0) && (pw0==0) && ( (pw2 ! =9) I I (pw6 ! =9) I I (pw3 ! =9) ) ) I I 
( (wetf lag==0) && ( (tempold [0] <temp_dew) && (tempold [0] >0 . 0) ) ) ) { 
wetf lag=l; 
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rhsum+= (float) rh ; 
rhinit= (float) rh ; 
twinit=tempold[0]  ; 

printfC'Wet  from  %d:00  7.d/7#d/7»d  with  T=  7#lf  ",hour, 
month, day , year ,tempnew[0] ) ; 
timeorg= (day*24+days [month] *24+hour) ; 

> 

if ( (wetf lag==l) && (pw0==0) && ( (pw2==9)&& (pw6==9) && (pw3==9) ) && 

(tempold [0] >=temp_dew) ) { 
wetflag=0; 

printf  ("until  7.d:00  7»d/7#d/7.d  ",  hour , month,  day , year) ; 

timef in= (day*24+days [month] *24+hour) ; 

dur at i on=t imef in-t ime o r g ; 

durcum+=duration ; 

numcum+=l ; 

printf  ("of  duration  */,d  hour  s \n ",  duration)  ; 
fprintf  (wetf ile, M,/»ld  */»d  7»f  7#f  \n" , time org, duration, 
(twinit+tempold[0])/2. ,rhinit) ; 

> 

> 

printf ("\n") ; 

printf  ("total  cumulative  time  of  wetness  is  7»ld  in  7.d  separate  events\n", 
durcum , numcum) ; 
rhave=rhsum/ (float )numcum ; 

printf  ("average  rh  before  start  of  rainfall  is  7*f  \n",rhave); 
f close (infile) ; 
f close (wetf ile) ; 
f close (frf ile) ; 
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